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ABSTRACT 


In  the  accelerating  flow  of  a  lighter  continuous  phase 
through  a  heavier  one,  small  non-uniformities  grow  into  large 
ones  due  to  the  Rayleigh-Taylor  Instability.  An  experiment 
exemplifying  the  large  bubble1  formation  due  to  Rayleigh-Taylor 
Instability  has  been  performed  and  simulated  using  the  PHOENICS 
8i  computer  code.  The  same  numerical  procedure  was  applied  to 
the  two-phase  flow  in  a  gun  barrel.  It  showed  that  the 
acceleration  provided  by  the  movement  of  the  projectile  can  cause 
initial  non-uniformities  to  grow  with  time. 
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In  the  accelerating  flow  of  a  lighter  continuous  phase 
through  a  heavier  one,  non-uniformities  grow  in  size  due  to  the 


Rayleigh-Taylor  Instability  (RTI)  forming  large  low-particle- 
density  'bubbles'  moving  through  high-particle-density  fluids. 

The  RTI  (Ref.  1)  has  for  a  long  time  been  recognised  when 
the  two  phases  are  liquids  with  differing  densities.  Several 
authors  have  demonstrated  the  growth  of  small  perturbances  into 
large  ones  in  experiments  where  for  example:  a  heavy  fluid  falls 
through  a  lighter  one  under  gravity  (Refs.  2,  3  and  4);  or  when  a 
lighter  fluid  rests  on  top  of  a  heavier  one  and  the  binary  fluid 
system  is  accelerated  downwards  at  a  rate  greater  than  that  of 
gravity  ( Ref.  5  and  6 )  . 

In  a  similar  manner,  a  gas  accelerating  through  a  bed  of 
particles  is  expected  to  act  on  non-uniformities  in  the  particle 
packing  distribution  causing  them  to  grow  in  size.  Such  a 
situation  is  encountered  in  guns  where  the  acceleration  is 
provided  by  the  projectile  movement  which  drags  the  gas  (and 
combustion  products)  with  it.  The  RTI  then  acts  on  non¬ 
uniformities  present  in  regions  of  the  gun  barrel  where  unburned 
or  partly  burnt  particles  exist  causing  them  to  grow. 

The  presence  of  bubbles  in  guns  would  cause  irregularities 
in  the  rate  of  combustion  of  the  particles.  The  transfer  of  heat 
in  the  flow  would  also  be  affected  since  the  mechanism  of  heat 
transfer  within  the  bubble  is  predominantly  convective  rather 
than  conductive.  Also,  the  bubbles  would  appreciably  affect  the 
radial  velocities  of  the  unburnt  or  partly  burnt  particles  thus 
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(a)  Devise  a  simple  experiment  which  will  exemplify  the 

formation  of  low  particle  density  bubbles  which  move  through 
high-particle-density  fluids; 

(b)  Devise  a  computer  model  which  simulates  the  observed 

phenomena  with  sufficient  accuracy; 

(c)  Use  this  model  with  appropriate  changes  in  the  initial  and 

boundary  conditions  and  in  the  material  properties  to 

explore  the  likelihood  of  Rayleigh-Taylor  instability  in 
guns . 


OUTLINE 


REPORT 


Section  4  of  this  report  describes  the  experimental  work 


carried  out  while  section  5  gives  the  results  of  the  simulation 
of  the  experiment  using  the  PHOENICS  84  computer  code.  In 
section  6,  a  computer  model  is  presented  for  simulating 


combustion  in  a  gun.  The  results  are  discussed  and 
recommendations  for  future  work  are  given  in  section  7. 


4 .  EXPERIMENTAL  WORK 


Experiments  are  considered  essential  since  they:  (l)  ensure 


that  all  the  processes  relevant  to  the  problem  in  question  are 
accounted  for  in  the  numerical  model;  and  (ii)  they  yield  basic 


experimental  data  for  comparison  with  the  numerical  predictions 


thus  providing  an  important  feedback  to  the  model  and  a  measure 
of  its  accuracy. 

In  the  present  work,  simple  experiments  were  performed  to 
investigate  the  Rayleigh-Taylor  phenomenon.  In  these 
experiments,  water  was  used  to  accelerate  through  a  bed  of  sand 
in  a  2D  glass  channel.  Using  water  is  not  unrealistic  since  in 
guns,  the  density  of  the  gas  producing  the  acceleration  is 
comparable  to  the  density  of  the  propellant  particles.  In  these 
experiments,  the  water  enters  the  bed  of  sand  through  a  porous 
side  wall.  This  to  some  extent,  simulates  in  the  laboratory,  the 
situation  in  guns  where  the  gas  producing  the  acceleration  comes 
from  a  distributed  source  which  is  the  burning  of  the  propellant 
particles . 

Two  different  experiments  have  been  performed.  In  the 
first,  the  interface  of  the  bed  of  sand  is  flat  while  in  the 
second  it  is  tilted.  In  the  first  experiment,  the  R-T  acts  upon 
inherent  non-uniformities  in  the  particle  packing  distribution 
while  in  the  second  an  additional  initial  disturbance  is  imposed 
through  tilting  the  interface.  These  two  experiments  are 
described  in  the  following  sections. 
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1.5cm. 


A  head 


tank  placed  approximately  3m  above  the  sand 


provides  the  pressure  head  required  to  push  the  water  through  the 
sand.  The  sand  particles  are  about  100pm  and  have  a  density  of 
2,500kg/m3.  Glass  beads  are  introduced  as  shwon  in  Figure  1  to 
break  down  the  turbulence  eddies  at  entry  and  to  produce  a  more 
uniform  flow  across  the  channel. 

A  video  system  capable  of  25  frames / second  was  used  to 
record  the  observations.  This  permitted  an  immediate  real  time 
or  frame-by-frame  replay  of  each  run.  The  photographs  presented 
in  this  report  were  taken  with  a  35mm  camera  fitted  with  a  motor 
drive  capable  of  5  frames / second . 

The  water  in  the  tank  was  dyed  blue  using  nigrosene.  The 
top  surface  of  the  water  was  therefore  always  discernible 
enabling  its  velocity  to  be  measured  at  all  times. 

4.1.2  Procedure 

Sand  is  poured  through  the  open  top  of  the  channel  to  the 
desired  height.  The  valve  (A)  (see  Figure  1)  is  then  slowly 
opened  to  allow  water  into  the  channel,  setting  the  particles 
into  motion  and  expelling  the  air  bubbles  trapped  within  the  sand 
bed.  When  the  water  reaches  the  top  of  the  channel,  valve  (A)  is 
closed  and  the  sand  allowed  to  settle.  Then  valve  (B)  is  opened 
to  drain  the  water  from  the  channel  but  not  enough  for  the  water 
level  to  fall  below  the  height  of  the  sand  bed.  This  operation 
is  repeated  several  times  until  all  the  air  traps  have  been 
expelled  from  the  system. 

The  settling  of  the  sand  after  the  water  inlet  valve  had 


been  closed.  never  left  a  flat  sand  interface. 


The  surface  of 


Eft 


m 


Fr.^r. 


the  sand  bed  sometimes  formed  a  sinusoidal  wave  whose  wavelength 
was  similar  to  the  width  of  the  channel  but  it  was  often  a  series 
of  random  small  waves  forming  an  irregular  wavy  sand  bed  surface. 
To  produce  a  flat  interface  whenever  recording  the  observations, 
the  side  of  the  channel  was  gently  tapped  and  although  these 
irregular  waves  at  the  top  of  the  sand  bed  disappeared  this  does 
not  mean  that  the  sand  'packing'  in  the  bulk  of  the  sand  bed  was 
uniform.  It  is  these  non-uniformities  that  get  the  Rayleigh- 
Taylor  Instability  started. 

4.1.3  Observations 

The  observations  made  are  shown  in  Figure  2.  At  t  =  os,  the 
initial  water  and  sand  levels  can  be  seen  prior  to  valve  opening. 
This  photograph  also  shows  in  the  background,  some  of  the 
equipment  that  was  used  in  the  experiment.  This  should  be 
ignored.  Upon  opening  the  valve,  the  bed  expands  particulately 
without  the  appearance  of  any  bubbles  in  it.  The  top  surface  of 
the  bed  curves  upwards  to  form  a  convex  surface  at  t  =  0.2s  but 
then  oscillates  and  at  t  =  0.6s  forms  a  concave  surface.  At  t- 
0.8s,  a  bubble  is  observed  in  the  main  body  of  the  bed.  This 
bubble  grows  with  time.  With  the  velocity  in  the  bubble  being 
higher  than  the  surrounding  mixture,  a  shear  instability  of  the 
Kelvin- Helmholtz  type  develops.  This  occurs  at  t  =  1.0s  with  the 
typical  mushroom  shapes  appearing.  Two  velocity  recirculation 
regions  or  vortices  form  near  each  of  the  walls.  These  regions 
grow  with  time  until  the  two  emanating  from  the  opposite  walls 


meet  in  the  middle  of  the  channel  at  t 


1.2s. 


Meanwhile  the 


leading  edge  of  the  bubble  moves  further  up  the  channel.  At  t  =  I 

1.4s,  a  somewhat  chaotic  behaviour  is  observed  where  the  sand 
which  has  been  brought  into  the  centre  of  the  channel,  is  being 

penetrated  by  the  water.  But  at  t  =  1.6s,  a  more  orderly  state  j 

is  restored  where  three  wavelengths  of  the  Kelvin-Helmholtz 
instability  can  be  seen  along  the  right  hand  wall. 

4 . 2  Experiment  with  a  Tilted  Interface 

The  Rayleigh- Taylor  instability  acts  upon  non-uniformities 
in  the  flow  and  causes  them  to  grow.  In  the  experiment  described 
in  section  4.1,  such  non-uniformities  are  present  but  they  are 
too  small  to  quantify.  These  non-uniformities  arise  due  to 
perhaps  slight  unintentional  deviations  of  the  channel  from  the 
vertical  or  from  the  uneven  rate  of  sedimentation  of  the  sand 
across  the  channel  width  (see  Section  4.1.2). 

This  however,  presents  a  problem  when  the  R-T  process  is  to 
be  modelled.  Unless  non-uniformities  are  specified  as  initial 
conditions,  the  calculations  will  not  show  any  growth  of  the  R-T 
instability.  An  experiment  was  therefore  carried  out  in  which  an 
initial  disturbance  was  imposed.  This  disturbance  was  introduced 

by  having  an  initially  tilted  sand  interface.  This  could  be 
easily  specified  as  an  initial  condition  in  the  numerical 
calculations  . 

4.2.1  Apparatus 

The  same  apparatus  used  for  the  flat  interface  experiments 
was  used  here.  However,  the  procedure  was  changed  to  produce  an 


initially  tilted  sand  interface. 


This  is  described  in  the  next 


sub- section . 


i>. 2. 2  Procedure 

Sand  is  poured  through  the  open  top  of  the  channel.  The 

o 

channel  is  then  tilted  by  approximately  10  to  the  vertical. 

Valve  (A)  (see  Figure  1)  is  then  slowly  opened,  letting  water 

into  the  test  section  and  expelling  the  air  bubbles  which  are 
trapped  in  the  sand  bed.  When  the  water  reaches  the  top  of  the 
test  section,  the  valve  (A)  is  closed  and  the  sand  allowed  to 
settle.  The  valve  (B)  is  then  opened  to  drain  the  water  from  the 

channel  but  not  enough  for  the  water  level  to  fall  below  the 

level  of  the  bed  of  sand.  This  procedure  is  repeated  several 
times  until  all  the  trapped  air  is  expelled  from  the  system.  The 
channel  is  then  set  vertical  but  the  interface  of  the  sand  bed 
remains  at  approximately  10°  tilt. 

As  in  the  flat  interface  experiment,  the  settling  of  the 
sand  was  not  uniform.  Ripples  always  appeared  at  the  surface  of 
the  sand  bed  and  these  were  eliminated  from  the  surface  of  the 
bed  by  gently  tapping  the  side  of  the  channel. 

4.2.3  Observations 

The  observations  made  are  shown  in  Figure  3.  At  t  =  Os,  the 
initial  tilted  sand  surface  can  be  seen.  The  equipment  showing  in 
the  background  should  be  ignored.  Upon  opening  the  valve,  the 
shallower  part  of  the  sand  bed  accelerates  at  a  faster  rate  (t  = 
0.2s!  and  a  bulge'  appears  in  the  left  hand  side  of  the  channel 
(t  =  0.4  and  0.6s).  At  t  =  1.0s,  the  first  traces  of  a  bubble 

1  1 


appear  which  grows  with  time.  Also,  recirculation  regions  form 
near  the  wall  due  to  K-H  instability.  The  K-H  instabilities  grow 
as  the  leading  edge  of  the  bubble  moves  further  up  the  channel  (t 


=  1.2s).  At  t  =  1.4s,  the  recirculation  zones  near  the  wall  are 
seen  to  be  stretched  and  at  t  =  1.6s,  a  'roll-over'  wave  on  the 
right-hand  side  wall  of  channel  is  observed. 


4 . 3  Discussion 

Two  experiments  have  been  described  above.  Despite  the 
difference  in  the  initial  conditions,  the  size  of  the  produced 


bubble  and  its  velocity  are  in  the  two  experiments  similar.  This 
can  be  explained  by  referring  to  the  basic  theory  of  the 
Rayleigh-Taylor  instability.  (Ref.  1).  Under  R-T  instability, 
the  growth  rate  of  a  wave  is  dependent  on  the  density  ratio  of 
the  two  phases,  the  acceleration  and  the  wavenumber  of  the 
disturbance  i.e. 


n  =  At  g  k 


where 


n  is  the  growth  rate 


At  is  the  Atwood  Number  given  by 

p  -  p 
2  1 

At  =  — - L 

n  +  n 
v  2  1 


and  k  is  the  wavenumber  given  by 


where  A  is  the  wavelength  of  the  disturbance. 

However,  in  any  system  certain  wavelengths  are  more  unstable 
than  others  and  therefore  grow  faster.  These  critical,  fast 
growing  wavelengths  are  determined  by  the  characteristics  of  the 
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system  such  as  viscosity  or  surface  tension.  In  the  present 

experiments,  it  is  the  viscosity  that  plays  this  important  role 

I 

I 

and  the  most  unstable  wavelengths  grow  faster  to  engulf  the 
smaller  ones.  In  the  end,  a  dominant  size  of  bubble  emerges 
which  is  restricted  by  the  width  of  the  channel. 

The  velocity  of  a  bubble  moving  up  a  channel  is  determined 
by  the  bubble  size  and  the  fluid  acceleration  (Refs.  5  and  7). 
Since  in  the  two  experiments,  the  bubble  size  is  similar,  its 
velocity  is  also  expected  to  be  similar. 

Another  feature  which  is  common  to  the  two  experiments  is 
the  appearance  of  the  velocity  recirculation  regions  near  the 
walls  due  to  the  K-H  instability.  This  arises  due  to  the 
velocity  difference  between  the  liquid  in  the  bubble  and  the  sand 
layers  on  its  borders.  In  the  two  experiments,  the  velocity  of 
the  bubble  is  similar. 

The  experiment  with  a  tilted  interface  is  suitable  for 
testing  the  predictions  of  PHOENICS  with  regards  to  the  growth  of 
the  R-T  instability.  These  predictions  are  given  in  the  next 
section . 

5 .  The  Numerical  Prediction  of  the  Experiment 

Despite  the  apparent  simplicity  of  the  experiment,  the 

numerical  task  is  quite  complex.  Not  only  is  the  situation  two- 

phase  and  two-dimensional,  it  also  involves  specifying  interphase 
transport  relationships.  Due  to  t£>e  physical  uncertainty 

regarding  these  processes  and  to  the  consequent  lack  of 

generalised  models,  a  numerical  study  was  carried  out  to 
determine  the  influence  and  relative  importance  of  viscosity,  the 
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interphase  friction  coefficient  and  the  effect  of  the  wall.  This 
investigation  is  reported  in  section  5.3.  In  the  following 
section,  the  two-phase  two-dimensional  balance  equations  solved 
using  the  PHOENICS  computer  code,  are  given. 


5 . 1  Governing  Equations 

In  the  absence  of  the  transfer  of  mass,  the  phase  mass 


conservation  or  continuity  equation  is  given  by:- 


-r—  p .  R  +  div(p  R  <  V  .  >  )  =  0 
Ot  1  1  111 


(5.1) 


where 

p  is  the  density 
R  is  the  volume  fraction 
<v>  is  the  velocity  vector 

and  the  subscript  i  refer  to  the  phase  in  question  i.e.  liquid  or 
solid . 


The  momentum  equation  is  given  by: 


aT  (p  R  ♦)  +  div  {R.  (p.  < V . >  ♦  -  p  grad  ♦)} 
o  t  1  1  1  11 


=  R.  (F  -  grad  p)  +  F.  ♦  F 
l  g  l  w 


(5.2) 


where 


t 

* 


2 


is  the  interphase  friction 

The  expressions  used  for  F.,  p  and  F  are  given  below  in  sections 

1  w 

5.3.1.  5.3.2  and  5.3.3  respectively.  First,  however,  a  brief 
descritpion  of  the  numerical  solution  procedure  embodied  in  the 
PHOENICS  81  code  is  given. 

5 . 2  THE  NUMERICAL  SOLUTION  PROCEDURE 

The  above  equations  were  solved  using  the  PHOENICS  81 
computer  code  (Refs.  8,  9  and  10). 

A  conventional  staggered  grid  is  used  where  the  velocities 
are  stored  at  the  centre  of  the  cell  faces  to  which  they  are 
normal  while  all  other  variables  are  stored  at  the  centres  of  the 
cells  themselves. 

The  velocity  locations  have  their  own  surrounding  cells, 
which  act  as  control  volumes  over  which  the  differential  momentum 
equations  are  integrated,  to  yield  the  corresponding  finite 
domain  equations  for  velocity.  Equations  for  the  other  dependent 
variables  employ  control  volumes  around  the  grid  points. 

The  result  of  the  integration  is  a  set  of  finite-domain 
equations  which  include  contributions  from  the  transient, 
convection,  diffusion  and  source  terms.  Upwind  differencing  is 
used  in  evaluating  the  convection  terms.  A  fully  implicit 
formulation  is  used. 

The  finite-domain  equations  are  solved  using  the  SIMPLEST 
and  IPSA  algorithms  (Ref.  10).  The  integration  proceeds  along 
the  Z  direction  (see  Figure  5)  from  the  bottom  to  the  top,  and  is 


repeated  until  convergence  is  achieved. 


Further  details  may  be 
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FIP  =  CFIPS  x  gi  x  R^  x  R^ 

where  the  subscripts  1  and  2  refer  to  the  light  and  heavy  phases 

respectively  and  CFIPS  is  the  interphase  momentum  transfer 

coefficient.  In  the  following,  CFIPS  took  the  values  1 0 6 ,  103 

2 

and  10  .  The  results  are  given  in  Figure  6.  These  show  the 

C 

volume  fraction  contour  plots  for  (a)  CFIPS  =  10  ,  (b)  CFIPS  = 

3  2 

10  and  (c)  CFIPS  =  10  at  0.2s  intervals  from  0  to  1.6s.  At  t  = 

0s,  the  bed  is  stationary.  Note  that  the  initial  disturbance  is 

hardly  discernible  at  the  surface  of  the  sand  bed. 

As  the  water  flows  in,  the  bed  expands  becoming  dilute  very 

quickly  in  the  lower  regions  of  the  domain.  In  all  three  cases  a 

bubble  forms  and  gradually  grows  with  time. 

6  3 

Changing  the  CFIPS  values  from  10  to  10  did  not  have  much 
influence  on  the  results.  However  for  CFIPS  =  1 0 2 ,  a 

considerable  difference  is  noted  with  the  heavier  phase  moving 
slower  up  the  channel  due  to  a  reduced  momentum  transfer  from  the 
lighter  phase.  Also  a  part  of  the  heavier  phase  remains 
undisturbed  in  the  lower  region  of  the  calculation  domain. 

5.3.2  The  Effect  of  Viscosity 

The  relationship  between  the  shear  stress  and  the  rate  of 
strain  for  two  phase  mixtures  is  not  linear.  It  depends  on  the 
heavier  phase  volume  fraction,  R^.  In  the  present  calculation 
the  following  relationship  (Ref.  12)  was  used:- 

M  =  Mf  exp  [k  R?]  ( 5 . A  ) 

where 
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k  =  2.5  ♦  U* 


d0’5 


where 

d  is  the  diameter  of  the  particles  in  microns 

and 

♦  is  the  shape  factor  which  is  equal  to  unity  for  spherical 
particles 

For  the  present  system,  k  =  3.9.  However,  changing  the  value  of 
k  was  found  to  appreciably  affect  the  results  of  the  calculation. 
In  the  following,  k  took  the  values  3.9,  10  and  15. 

The  results  are  given  in  Figure  7.  These  show  the  volume 
fraction  contour  plots  for  (a)  k  =  3.9  (b)  k  =  10  and  (c)  k  =  15, 
when  CF IPS  =  103. 

As  the  water  flows  in,  the  bed  expands  becoming  dilute  very 
quickly  in  the  lower  region  of  the  domain.  For  the  higher  values 
of  k,  a  more  pronounced  bubble  is  formed  with  a  velocity 
recirculation  region  appearing  near  the  wall  at  later  times. 

It  should  be  noted  that  the  results  shown  here  are  not  grid 
independent.  A  further  increase  in  the  number  of  cells 
especially  in  the  radial  direction  will  affect  the  results. 
Improving  the  accuracy  of  the  results  is  however,  not  thought  to 
be  worthwhile  in  this  part  of  the  study. 

5.3.3  The  Effect  of  the  Wall 

The  wall  shear  stress  is  given  by:- 

t 

w 


p  grad  V 


Since  the  viscosity,  g  features  in  the  wall  shear  stress 


expression,  it  was  decided  to  study  the  isolated  effect  of  the 
wall.  A  wall  coefficient,  C,  was  introduced  to  give:- 

(5.5) 

Keeping  g  constant,  C  was  given  the  values  1  and  10. 

The  results  are  given  in  Figure  8.  These  show  the  volume 


t  =  Cg  grad  V 
w 


fraction  contour  plots  for  (a)  C  =  1  and  k  =  3.9,  (b)  C  =  10  and 

.3 


k  =  3.9  and  (c)  =  10  and  k  =  10  when  CFIPS  =  10" 


As  the 

water 

flows 

in 

.  the  bed 

expand 

s  becoming 

dilute  quite 

quickly.  A 

bubble  forms 

quicker 

when  C 

=  10  but  at 

later  times 

it  can  be 

seen 

that 

th 

e  effect 

of  the 

wall  is 

only  local. 

Increasing 

k  to 

10 

can 

be  seen 

not  t 

o  affect 

the  results 

appreciably . 

5.3.4  Conclusions 

Preliminary  calculations  have  established  that  unless  an 
initial  disturbance  is  specified  in  the  calculations,  a  bubble 
will  not  form.  With  a  small  initial  disturbance,  a  parametric 
study  has  been  conducted. 

The  results  from  the  study  have  shown  that  interphase 
friction  and  the  mixture'  viscosity  specification  are  important 
to  the  formation  and  behaviour  of  the  bubble.  The  effect  of  the 
wall  was  seen  to  be  local. 


5 . 4  The  Prediction  of  the  Experiment  with  a  Tilted  Interface 

The  experiment  with  a  tilted  interface  has  been  simulated  as 
a  two-dimensional  two-phase  problem.  The  continuity  and  momentum 


conservation  equations  for  each  of  the  phases  (see  Section  5.1) 


nave  been  solved  using  PHOtNlCS.  The  interphase  friction 
description  given  in  equation  5.3  was  used  with  CFIPS  =  103.  The 
viscosity  as  given  in  equation  5.4  was  used  with  k  =  10.  The 


effect  of  the  front,  back  and  two  side  walls  were  taken  into 
account  with  the  constant  C  in  equation  5.5  set  equal  to  1. 

The  rate  of  expansion  of  the  part-fixed  moving  grid 


corresponded  to  the  mass  flow  rate  into  the  channel  such  that  the 
water  level  in  the  experiment  coincided  with  the  moving  boundary 
of  the  grid  at  all  times.  The  grid  consisted  of  30  cells  in  the 


lateral  direction  and  60  cells  in  the  axial  direction.  The  time 
step  used  was  0.005  seconds. 

The  experimental  observations  alongside  the  volume  fraction 


( R 2 )  contour  plots  and  the  light  and  heavy  phase  velocities  (l^ 
and  W.,)  are  shown  in  Figure  9.  At  t  =  os,  the  initially  tilted 
interface  can  be  clearly  seen.  As  the  water  enters  the  domain 


from  the  back,  the  shallower  left-hand  side  of  the  bed  is 
accelerated  at  a  faster  rate.  At  t  =  0.4  and  0.6s,  the 
experiment  shows  a  bigger  bulge'  believed  to  be  due  to  three- 
dimensional  effects  in  the  plane  perpendicular  to  the  paper.  At 
t  =  0.8s,  the  bubble  appearing  in  the  experiments  is  seen  in  the 
channel.  The  bubble  continues  to  grow  and  velocity  recirculation 


regions  due  to  a  Kelvin-Helmholtz  instability  can  be  seen  to  have 
started  to  form  at  t  =  1.0s.  These  continue  to  grow  at  t  =  1.2s 
while  the  leading  edge  of  the  bubble  moves  further  up  the 


channel.  The  mushroom  shape  appearing  in  the  experiment  at  t  = 
1.2s  appears  in  the  calculations  at  t  -  1.4s.  This  continues  to 


move  up  the  channel  at  t  =  1.6s. 


The  model  has  predicted  in  the  experiment  the  initial  uneven 


acceleration  of  the  different  parts  of  the  bed,  the  formation  and 
growth  of  the  bubble  and  the  velocity  recirculation  regions  that 
appear.  It  has  however  overestimated  the  rate  at  which  the  sand 
moves  up  the  channel  as  a  whole.  This  is  believed  to  be  due  to 
the  somewhat  simple  relationships  used  for  the  effective 
viscosity  and  the  interphase  friction. 

It  should  be  emphasised  that  the  results  obtained  here  have 
not  been  tested  for  grid  independence.  In  the  author's 

experience  with  this  calcualtion,  the  refinement  of  the  grid  is 
important  especially  in  the  radial  direction. 

Finally,  it  should  be  said  that  there  is  a  possibility  that 
an  oversimplification  has  been  committed  in  treating  the  problem 
as  two-dimensional.  Although  a  three-dimensional  calcualtion  can 
be  carried  out,  the  effort  involved  cannot  be  justified  for  the 
relatively  small  improvement  that  this  would  make. 

5 . 6  CONCLUSIONS 

From  this  part  of  the  work,  the  following  conclusions  can  be 
drawn : - 

1.  The  experiment  described  here  has  shown  that  when  a  lighter 

continuous  phase  accelerates  through  a  heavier  one,  small 
non-uniformities  grow  into  large  ones  due  to  the  Rayleigh- 

Taylor  Instability. 

2.  Despite  the  absence  of  combustion  from  the  experiment,  it 


simulates  the  hydrodynamics  of  the  gun-barrel  situation 
specifically  with  regard  to  the  acceleration  of  a  light 
fluid  through  a  heavier  one. 

3.  The  two-phase  numerical  model  employed  has  predicted  fairly 

well  the  Rayleigh-Taylor  Instability.  The  discrepancies 
between  the  predictions  and  the  experimental  observations 
are  attributed  to  uncertainties  about  the  inter-phase 

friction  and  effective- viscosity  laws  which  have  been  used. 

4.  The  quality  of  the  agreement  between  predictions  and  the 

observations  gives  credibility  to  the  predictions  which  the 

same  computational  procedure  produces  when  it  is  applied  to 
a  gun-barrel  simulation  involving  combustion. 

The  remaining  part  of  this  report  dedicates  itself  to  the 

prediction  of  the  gun-barrel  situation  especially  with  regard  to 
the  growth  of  large  non-uniformities  from  smaller  ones. 


It  has  been  demonstrated  above  that  the  acceleration  of  a 

fluid  containing  a  dispersed  heavier  phase  leads  to  large  bubble 

formation.  a  consequence  of  the  Rayleigh-Taylor  Instability.  In 
guns,  initial  particle  packing  non-uniformities  or  those  that 

develop  as  the  ignitor  discharges  the  hot  gas  into  the  propellant 
bed  may  prove  to  be  the  source  of  volume  fraction  discontinuities 
resulting  from  the  Rayleigh-Taylor  instability.  The  objective  in 


this  part  of  the  work  is  to  establish  whether  low-particle- 


density  regions  or  'bubbles'  can  occur  under  conditions  which  are 
true  of  guns. 

If  R-T  bubbles  do  occur  in  guns,  then  irregularities  in  the 
rate  of  combustion  of  the  propellant  particles  can  be  expected. 
Also,  the  bubbles  would  appreciably  affect  the  radial  velocities 
of  the  unburnt  or  partly  burned  particles  thus  influencing  the 
heat  transfer  to  and  erosion  of  the  barrel  walls. 

In  what  follows,  two  calculations  are  presented.  The  first 
simulates  a  gun-like  situation  where  particles  contained  in  a 
cylindrical  domain  ignite,  which  raises  the  pressure  within  the 
domain  causing  it  to  expand.  The  expanding  domain  simulates  the 
movement  of  the  projectile  in  guns.  The  results  from  such  a 
calculation  are  described  and  discussed. 

The  second  calculation  is  identical  to  the  first  except  that 
a  small  non-uniformity  in  the  volume  fraction  distribution  is 
introduced.  The  results  from  this  calculation  are  described  and 
compared  with  the  previous  calculation  especially  with  regard  to 
the  effect  that  the  initial  non-uniformity  has  had. 

It  is  worth  noting  that  the  calculations  presented  here  do 
not  provide  quantitative  results  of  the  gun  barrel  problem.  They 
are  aimed  more  at  demonstrating  that,  with  the  acceleration 
provided  by  the  movement  of  the  projectile,  the  non-uniformities 


that  may  be  present  grow  in  size. 


The  basic  features  of  the  gun  barrel  are  outlined  in  Figure 
10.  A  cylindrical  domain  is  considered,  enclosed  by  the  gun 
barrel  and  the  base  of  the  projectile,  and  containing  a  solid 
propellant  and  a  gas.  The  -olid  propellant  is  assumed  to  be 
spherical  particles. 

Ignition  is  provided  by  the  inflow  of  hot  gases  at  the  base 
of  the  domain.  The  gas  is  forced  into  the  propellant  bed  which 
causes  a  compaction  of  the  granular  bed  near  the  entrance  region 
and  also  heats  up  the  nearby  granular  propellants  to  ignition. 
The  ignited  propellants  give  off  more  hot  gases  which  are  pushed 
forward  by  the  pressure  gradient  to  ignite  more  propellants. 
Thus  a  pressure  gradient  is  created  inside  the  combustion  chamber 
and  the  accelerated  gaseous  products  cause  the  projectile  to 
move . 


The  modelling  of  such  unsteady  two-phase  flow  phenomena 
requires  the  solution  of  the  unsteady  two-phase  gas  dynamic 
equtaions  as  well  as  the  utilisation  of  associated  empirical 
correlations  for  interphase  friction  and  heat  transfer 
coefficients  and  burning  rate  laws  for  the  rate  of  interphase 
mass  transfer. 

6 . 3  The  Dependent  and  Independent  Variables 

The  dependent  variables  of  the  problems  are:  the  velocities 
of  the  gas  and  particles  in  the  axial  and  radial  directions,  W 
W^,  V ^  and  :  gas  and  particle  volume  fractions  and  R  ;  and 
enthalpies  of  gas  and  particles  h^  and  h  . 
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The  independent  variables  are:  the  axial  and  radial 

distances  z  and  r;  and  the  time  t. 

6 . 4  The  Partial  Differential  Equations 

The  following  set  of  governing  equations  describe  the  change 
of  mass,  momentum  and  energy  for  each  of  the  gas  and  solid 
phases . 


6.4.1  The  Mass  Conservation  Equation 
Gas  Phase 
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Particle  Phase 
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P1  and  p^  are  the  densities  of  the  gas  and  particle 
respectively 

Vj ,  V  ,  W1  and  are  the  velocities  of  the  gas  and  particle 
in  the  radial  and  axial  directions 


and 

m"  is  the  rate  of  mass  transfer  per  unit  volume  from  the 
solid  phase  to  the  gaseous  phase  due  to  gasification 
of  the  soli>‘  particles 

The  volumetric  fractions  are  related  by  the  space-sharing 


equation : - 


6.4.2  The  Conservation  of  Momentum  Equations 

The  momentum  equations  for  a  transient  two-phase  flow  are:- 


3  1  3  3 

—  (p  R.  ♦)  +  -  (r  R  p.  V.  ♦)  ♦  —  (R.  p.  W.  4) 
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where 

4  stands  for  V  ,  V2,  Wi  and  W2; 

r  is  the  diffusion  coefficient 

represents  source/sink  of  momentum  due  to  pressure 
gradients,  gravitational  forces  etc... 
and  the  subscript  i  refers  to  the  phase  in  question  i.e.  gaseous 
or  soliC. 

For  the  application  of  the  model  considered  in  this  part  of 
the  work,  diffusion  effects  are  considered  negligible  =  0). 

The  source  terms  on  the  RHS  of  equation  6.4  are  given  in 
Table  1.  The  effect  of  the  wall  has  been  neglected.  Note  that 
the  pressure  gradient  term  is  written  as  R  grad  p  as  opposed  to 
grad  (p  R).  The  form  used  in  the  present  work  is  the  correct  one 
(Ref.  13). 

In  the  present  formulation,  the  two  phases  share  the  same 
pressure.  However,  an  additional  intergranular  force'  term, 
grad  (R2  t),  appears  in  the  solid  phase  equation.  This  describes 
the  extra  stress,  t,  sustained  by  the  solid  phase  as  its  volume 
fraction  approaches  the  physically  attainable  limit.  This 
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intergranular  stress  is  a  function  of  the  volume  fraction 
The  m21  <V^>  terms  describe  the  rate  of  change  of 
due  to  the  motion  of  the  gasifying  particles. 


momentum 


6.4.3  The  Conservation  of  Energy  Equations 
Gas  Phase 
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Particle  Phase 
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where 


=  P  ♦  t 


q  is  the  rate  of  heat  transfer  from  the  gas  to  the 


particles 


=  h2  +  hc 


where 


hc  is  the  heat  of  combustion  of  the  solid  particles. 


In  the  above,  the  rate  of  heat  transfer  within  the  solid  phase  by 
conduction  has  been  neglected.  The  rate  of  heat  transfer  to  the 
gun  barrel  wall  has  also  been  neglected. 


Equation  6.1  to  6.6  form  a  system  of  nine  coupled,  non- 


linear,  partial  differential  equations.  These  have  been  solved 
using  the  PHOENICS  computer  code. 


6 . 5  Auxilliarv  Relations 

To  complete  the  mathematical  specification  of  the  problem 
constitutive  relations  are  needed.  Some  of  the  relations 
presented  below  are  simple  and  need  improvement  if  accurate 
results  are  to  be  obtained.  These  however  are  not  essential  to 
the  Rayleigh-Taylor  problem. 


6.5.1  Equation  of  State 

The  equation  of  state  for  a  perfect  gas  undergoing  an 
lsentropic  process  was  used.  This  is  given  by:- 

pp  a  =  const .  (6.7) 

A  more  representative  equation  of  state  would  be  the  Nobel- 
Abel  or  Clausius  equation  (Ref.  K)  which  takes  into  account  the 
molecular  volume. 

6.5.2  The  Intergranular  Stress 

A  granular  packed  bed  under  compressive  load  can  be  further 
compacted.  A  measure  of  this  compaction  is  the  volume  fraction 
of  the  heavy  phase,  R^.  There  is  however,  a  limiting  maximum 
compaction  depending  on  the  particle  shape,  properties  and  size 
distribution.  In  the  case  of  unisized,  incompressible  and 
spherical  solids,  the  maximum  compaction  corresponds  to 
approximately  R2  =  0.75.  If  the  particles  are  of  different  sizes 


in  the  bed,  higher  compaction  is  possible. 
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e  and  n  are  constants. 

The  production  rate  of  gases  from  the  solid  particles  is 
therefore  given  by:- 

m21  =  b  As  (6.10) 

where 

A  is  the  total  surface  area  of  particles  in  a  control  cell 
s 

calculated  for  spherical  particles  from 

6  R  Vol 

A  =  - - 

S  0 

o 

where 

Vol  is  the  control  cell  volume 
and 

0q  is  the  initial  diameter  of  unburnt  particles 

It  should  be  noted  that  the  present  calculations  do  not 
account  for  the  reduction  in  particle  size  due  to  combustion. 
Consequently,  the  rate  of  gas  generation  is  therefore 
overestimated . 

An  important  improvement  to  the  present  model  would  employ 
the  shadow'  method  for  particle  size  calculation  (Ref.  18). 


I 

The  particle  is  assumed  to  ignite  when  its  surface 

temperature  reaches  a  specified  value  i.e. 

b  =  o  T  <  T 

s  ignition 

n  (6.11) 

-) . ,,  , 

I  s  ignition 

atm  J 

The  particle  surface  temperature  is  determined  from  equation  6.15 
below. 

6.5.6  The  Interohase  Heat  Transfer  Coefficients 

Although  the  equations  solved  for  the  transport  of  heat 
between  the  gas  and  particle  phases  are  those  of  the  phase 
enthalpies  (equations  6.5  and  6.6),  it  is  convenient  to  think  in 
terms  of  the  temperatures  T^  and  T^  by  introducing  the  specific 
heats  Cj  and  for  the  gas  and  particle  phases  respectively. 
For  simplicity,  and  are  assumed  equal  (=  C). 


Central  to  the  following  treatment  is  the  concept  of  an 
interface  between  the  two  phases  with  temperature  T  Then 


q ,  -  a  ( T 

-  T  ) 

(6.12) 

Is  1  1 

s 

q  =  a„  ( T 

-  T,) 

(6.13) 

s2  2  s 

2 

where  the  subscripts  1 , 

2  and  s  refer  to  the  gas, 

solid  and  the 

interface  respectively. 

q  is  the  rate  of 

heat  transfer  and  a  and  a 

are  the  heat 

transfer  coefficients  multiplied  by  the  interface  area  through 
which  the  transfer  occurs. 

An  energy  balance  over  the  control  volume  enclosing  the 
interface  yields : 
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heat  coming  into  the  control  volume 


*  V  *  "2.  CT2 

heat  going  out  of  the  control  volume 


=  °S2  *  "2,  CTS 

generation  =  rn  h 
21  c 


where 


m  is  given  in  equations  6.10  and  6.11 


and 


h^  is  the  heat  of  combustion  of  the  solid  particles 
Therefore 


q  -  q 

s2  is 

1  m21 

Ch  - 
c 

C  ( T  - 
s 

V’ 

Combining 

equations  6. 

12.  6. 

13  and 

6 . H  yields 

T  = 
s 

Ca2T2 

+  aiT1 

*  "21 

( h  ♦ 
c 

ct2)]/a 

<1. 

=  31/A 

Ca2(T1 

-  T2> 

+  *21 

C(T1  -  T2 ’  - 

«s2 

1  a2M 

Cai(T1 

-  T2  ’ 

+  *21 

h  ] 
c 

where 

A  -  3  t 

+  3  2  + 

m21C 

(6  .14) 


(6.15) 


21c' 


The  heat  transfer  coefficient  per  unit  area, 


1 


is 


calculated  from  the  Oenton  (Ref.  19)  correlation  as  modified  by 
Eckert  and  Drake  (Ref.  20)  given  by:- 
D 


Nu  = 


1 


-  0,.  Re0'7  Pr°'3 


1 


where 


k, 

is 

the 

thermal 

conductivity 

of  the  gas 

0 

1  s 

the 

particle 

diameter 

Re 

1  s 

the 

Reynolds 

Number  given 

by :  - 

Re 

-  WZ)R10/(P1 

and 
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Pr  is  the  Prandtl  Number  given  by : - 
Pr  =  Mt  C/kt 


6.5.7  The  Movement  of  the  Pronectile 

The  projectile  velocity  and  acceleration  were  calculated  at 
every  time  step  from  the  computed  force,  F,  acting  on  its  base. 

F  =  (pf  '  Pb)A 

where 

P  =  projectile  frictional  pressure 

P  =  calculated  pressure  at  the  base  of  the  projectile 

b 

and 

A  :  projectile  cross-sectional  area 

6 . 6  The  Problem  Considered 

The  problem  considered  is  that  of  predicting  the  time- 
dependent  two-phase  two-dimensional  heat  transfer  and  combustion 
processes  occunng  in  a  gun  barrel  as  illustrated  in  Figure  10. 
Table  2  summarizes  the  input  data  used  and  the  initial  conditions 
employed . 

At  t  -  Os,  hot  gases  flow  into  the  domain.  The  heat 
transfer  to  the  particles  raises  their  surface  temperature  until 
the  ignition  temperature  is  reached  initiating  combustion. 
Further  hot  gases  are  produced  as  a  result  of  combustion  and  the 
pressure  builds  up  inside  the  domain.  The  temperature  also 
increases  igniting  more  particles.  When  the  pressure  is 

sufficiently  high,  the  projectile  starts  to  move. 


Two  calculations  are  presented  below. 


The  first  computes 


variations  with  time  along  the  axial  direction  only.  In  the 
second  calculation,  a  non-uniformity  in  the  volume  fraction 
distribution  is  introduced.  The  results  from  the  two 
calculations  are  compared  to  demonstrate  the  effect  of  the 
initial  non-uniformity  on  the  results. 


6 . 7  Computational  Details 

The  gun  barrel  shown  in  Figure  10  is  divided  into  a  number 
of  annular  cells  in  the  axial  and  radial  directions.  An  example 
is  given  in  Figure  11  which  shows  two  cells  in  the  axial  z 
direction  and  three  cells  in  the  radial  y  direction. 

A  numerical  test  was  performed  to  establish  the  dependence 
of  the  results  on  the  number  of  cells  in  the  axial  direction.  In 
a  10  calculation,  and  for  a  fixed  time  interval  of  0.04ms,  the 
number  of  ceils  was  varied  from  ten  to  thirty  keeping  all  other 
properties  unchanged.  The  dependence  of  the  pressure  on  the 
number  of  cells  is  shown  in  Figure  12.  It  can  be  seen  that 
little  accuracy  can  be  gained  from  increasing  the  number  of  axial 
cells  beyond  thirty.  In  the  following  calculations,  30  and  16 
cells  were  used  in  the  axial  and  radial  directions.  The  time 
interval  used  was  0.02ms. 

Following  the  projectile  movement,  the  grid  was  allowed  to 
expand  in  the  axial  direction  at  every  time  step  subject  to  the 
conditions  outlined  in  section  6.5.7  above.  The  lower  twenty 
cells  corresponding  to  the  first  0.13m  from  the  bottom  of  the 
barrel  were  not  allowed  to  move. 


Vc 


6 . 8  RESULTS 


6.8.1  Simulation  of  the  gun  situation 

The  results  from  are  shown  in  Figures  13  to  17.  The  axial 
distribution  histories  for  the  pressure,  gas  and  particle 
velocities,  particle  volume  fraction  and  temperature  are  given. 

In  figure  13,  the  axial  pressure  distribution  is  shown  at 
t=0.5,  1.5,  2.5,  3.5,  4.5,  5.5  and  6.5ms.  At  t=0.5ms,  the 
pressure  near  the  base  of  the  barrel  rises  sharply  due  to  the 
inflow  of  hot  gases  which  ignite  the  particles  in  this  region.  A 
steep  pressure  gradient  is  obtained  which  diminishes  with  time. 
The  pressure  builds  up  with  time  and  at  t=2.5ms  it  exceeds  1  M  Pa 
the  point  at  which  the  projectile  starts  to  move.  Despite  the 
movement  of  the  projectile  (at  t=2.16ms)  the  pressure  continues 
to  rise  due  to  particle  gasification  until  it  reaches  a  maximum 
of  around  20  M  Pa  at  t=6.5ms. 

Figures  14  and  15  show  the  gas  and  particle  axial  velocity 
distributions  at  the  same  times  as  those  for  the  pressure.  A 
peak  is  observed  in  both  the  gas  and  particle  profiles  up  to 
t=2.5ms  driven  by  the  pressure  gradient.  Beyond  t=2.5ms,  a  peak 
in  the  gas  velocity  profile  is  obtained  behind  the  pro3ectile  as 
it  drags  the  gas  with  it.  A  maximum  gas  velocity  of  about  300m/s 
is  obtained  at  t-6.5ms.  The  particle  velocity  is  much  lower  due 
to  the  higher  inertia  that  the  particles  possess.  A  maximum 
particle  velocity  of  about  130m/s  is  obtained  at  t=6.5ms. 

The  particle  volume  fraction  (R2)  distribution  history  is 
shown  in  figure  16.  Near  the  base  of  the  barrel,  the  R2  values 
decrease  quite  rapidaly  due  to  the  burning  of  the  particles  and 
their  forward  motion  as  they  are  carried  by  the  gas.  The  effect 
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of  the  latter  causes  a  compaction  of  the  particles  in  other 
regions  of  the  Darrel.  This  is  clearly  seen  in  figure  16  where 
the  R2  values  increase  beyond  the  initial  value  of  0.5.  This 
effect  is  however  exaggerated  here  due  to  the  neglect  of  the 
intergranular  stress  which  acts  to  resist  this  compaction  (see 
section  6.5.2).  Near  the  base  of  the  projectile,  the  sharp 
decrease  in  the  R2  values  for  t>2.5ms  is  caused  by  the  movement 
of  the  projectile. 

Figure  17  gives  the  gas  temperature,  T1,  axial  distribution 
history.  At  all  times  near  the  bottom  of  the  barrel,  T1  remains 
at  2000K  which  is  the  temperature  of  the  injected  gas  specified 
as  an  initial  condition.  This  temperature,  however,  decays 
rapidly  with  axial  distance  downstream.  With  time,  T1  builds  up 
gradually  in  the  top  part  of  the  barrel  due  to  heat  transfer  by 
convection  combined  with  combustion  -  generated  heat.  The 
particle  temperature  distribution  history  follows  closely  that 
for  T 1  due  to  the  high  heat  transfer  coefficient,  a^,  which  was 
used . 

Finally,  the  projectile  velocity  and  acceleration  are  shown 

in  Figure  18.  The  projectile  starts  moving  at  t=2.16ms  and 

2 

accelerates  rapidly  to  about  12000  m/ s  at  t=6  ms  after  which  the 
acceleration  decreases.  The  velocity  at  t=6ms  is  about  75m/s. 
In  6.5ms,  the  acceleration  of  the  projectile  produces  a  73. 5Z 
expansion  in  the  length  of  the  domain. 


6.8.2  The  effect  of  an  initial  non-umformitv 

An  initial  volume  fraction  non-uniformity  was  introduced  in 
the  calculation.  This  non-uniformity  consisted  of  a  cell  in 
which  the  particle  volume  fraction,  R2,  was  reduced  from  0.5  to 
0.1  whxle  its  two  neighbouring  cells  in  the  radial  direction  had 
R2  =  0.7.  This  disturbance  occupied  only  one  cell  in  the  axial  Z 
direction  at  Z/L  =  0.33  and  r/R  =  0.64.  Since  the  volume  of  the 
computational  cells  is  everywhere  constant.  the  total  amount  of 
propellant  mass  is  therefore  unchanged.  With  the  exception  of 
the  three  above  mentioned  cells,  the  R2  values  were  everywhere 
equal  to  0.5. 

The  effect  of  this  non-uniformity  on  the  pressure  and 
temperature  has  not  been  significant.  The  pressure  and 
temperature  distribution  histories  are  very  similar  to  figures  13 
and  17  with  no  significant  variations  in  the  radial  direction. 

The  gas  velocity  axial  distribution  history  given  in  figure 
19  shows  a  sudden  increase  in  the  gas  velocity  at  t  =  0.5ms  and 
Z/L  =  0.33.  This  coincides  in  position  with  the  initial  non- 
uniformity  and  is  typical  of  the  RTI.  For  the  RTI  accelerates 
unevenly  the  parts  of  the  flow  which  contain  volume  fraction  non- 
uniformities .  The  acceleration  for  t<2.5ms  is  not  provided  by 
the  movement  of  the  projectile  (as  it  has  not  yet  started  to 
move)  but  is  due  to  the  pressure  gradient  within  the  barrel  (see 
figure  13). 


The  uneven  axial  gas  velocity  distribution  is  evident 
throughout  the  calculation  although  figure  19  does  not  show  it. 
In  figure  20.  the  radial  gas  velocity  distribution  at  la)  Z/L  = 
0.33  and  fbi  Z/L  -  0.88  is  given.  At  Z/L  =  0.33  the  non-uniform 
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gas  velocity  is  seen  to  be  at  its  greatest  at  t  -  0.5ms  after 
which  it  diminishes  due  to  the  forward  movement  of  the  volume 
fraction  non-uniformity  and  the  reduction  in  its  magnitude  due  to 
the  effects  of  compaction.  At  Z/L  =  0.88  (figure  20  (b)),  the 
uneven  velocity  profiles  are  only  seen  at  later  times. 

The  volume  fraction  axial  distriDution  history  is  shown  in 
figure  21.  For  clarity  of  presentation,  the  R2  profiles  have 
been  given  in  figure  21  (a)  at  t=0.5,  1.5  ad  2.5ms  and  in  figure 
21  (b)  at  t-3.5,  4.5,  5.5  and  6.5ms.  At  t-0.5ms,  the  volume 
fraction  non- uniformity  is  seen  to  extend  from  Z/L  =0.35  to  Z/L  = 
0.45  thus  occupying  about  10Z  of  the  domain  length.  Its 
magnitude  however  ranges  between  R2  =  0.15  to  0.5.  At  subsequent 
times,  the  magnitude  of  the  non-uniformity  is  seen  to  diminish 
but  it  grows  in  size.  For  at  t=4.5ms  (figure  21  ( b ) )  the  non¬ 
uniformity  occupies  about  30Z  event  though  its  magnitude  ranges 
between  R2=0.7  and  0.8.  The  decrease  in  the  magnitude  of  the 
non-uniformity  is  due  to  the  effects  of  compaction  but  the 
increase  in  its  size  is  due  to  RTI. 

The  significance  of  the  RTI  can  be  appreciated  when 
comparing  the  volume  fraction  profiles  from  this  calculation  with 
the  corresponding  ones  in  the  absence  of  the  initial  non¬ 
uniformity.  This  is  done  in  figure  22  for  (a)  t=0.5ms,  (b) 
t=2.5ms,  (c)  t=4.5ms  and  (d)  t=6.5ms.  The  increase  in  the  size 
of  the  non-uniformity  can  be  clearly  seen  although  at  t=6.5ms  it 
is  affected  by  the  movement  of  the  projectile. 

The  volume  fraction  radial  profiles  at  Z/L  =0.33  and  Z/L 
=0.88  are  given  in  figure  23  (a)  and  (b).  The  non-uniformity  is 


seen  not  to  grow  in  the  radial  direction. 


6.8.3  Discussion 

The  calculation  presented  in  section  6.8.1  has  produced 
realistic  results  from  a  simplified  model  for  the  combustion  of 
propellant  particles  in  a  gun  barrel.  The  inflow  of  the  hot 
gases  causes  the  particles  to  ignite  raising  the  pressure  and 
temperature,  thus  igniting  more  particles  and  causing  the 
projectile  to  eventually  move.  The  inflow  of  hot  gases  pushes 
the  particles  away  from  the  injection  region  and  causes  them  to 
compact  everywhere  else.  The  neglect  of  intergranular  stress 
from  these  calculations  has  given  rise  to  high  compaction  and 
subsequently  high  volume  fraction  values. 

The  introduction  of  a  small  non-uniformity  has  not  affected 
the  pressure  nor  the  temperature  significantly.  But  the  initial 
volume  fraction  non-uniformity  was  seen  to  increase  in  size  as  it 
moves  along  the  barrel  despite  a  decrease  in  its  magnitude.  The 
magnitude  of  the  non-uniformity  is  believed  to  have  been  strongly 
affected  by  the  compaction  of  the  particles  which  has  been 
exaggerated  through  the  neglect  of  intergranular  stress.  But  the 
RTI  is  still  evident  even  though  it  is  superimposed  on  top  of 
other  effects.  A  sketch  of  the  author's  view  of  what  happens  is 
given  in  figure  24. 

i 

i 

i  6.9  Conclusions 

i 

I 

i  From  this  part  of  the  work,  the  following  conclusions  can  be 

drawn . 


1.  The  simulation  of  a  simplified  gun  barrel  has  produced 
realistic  results  using  a  numerical  procedure  which 
was  previously  shown  to  predict  the  RT I  satisfactorily. 

2.  A  calculation,  in  which  a  small  initial  non-uniformity 

in  the  particle  volume  fraction  has  been  introduced, 
demonstrated  that  these  do  grow  as  a  consequence  of  the 
Rayleigh-Taylor  Instability. 

3.  The  growth  of  the  non-uniformity  and  the  common  features 
that  this  calculation  has  with  the  gun  barrel  situation 
gives  credibility  to  the  notion  that  the  possibility  of  RT I 
in  guns  deserves  attention. 

T .  Future  Work 

The  model  used  for  simulating  the  combustion  within  the  gun 
barrel  can  be  improved  to  provide  realistic  predictions.  These 
improvements  would  involve  better  expressions  for  the  equation  of 
state,  the  interphase  friction  and  the  intergranular  stress.  A 
particle  size  calculation  should  also  be  introduced.  Recommended 
expressions  have  been  given  in  section  6.5. 
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TABLE  1 


n  the  momentum  eaua 


cm]<o  N  kO 


Parameter 


Value 


Physical  Properties 


Propellant  density  1500kg/m3 

Specific  heat  ratio  of  gas,  a  1.4 

Ignition  temperature  400K 

Chemical  energy  of  propellant,  h  4  MJ/kg 

Specific  heat  of  gas  and  particle,  C  2000  J/k^  K 

Mass  in  flow  rate  from  ignitor  40kg/s  m 

Temperature  of  gas  from  ignitor  2000K 


Constitutive  Relations 


Propellant  burning  rate  proportionality  2 

constant,  e  0.2kg/m  s 

Propellant  burning  rate  index,  n  0.9 

Interphase  friction  parameter,  C  (eqn.  6.9)  100 

Interphase  heat  transfer  parameter  on 

particle  side,  a?  20W/K 

Initial  Conditions 

5  2 

Pressure  1x10  N/m 

Bulk  temperature  of  solid  particles  294K 

Temperature  of  gas  294K 

Volume  fraction  of  solid,  0.5 

Velocities  of  gas  and  particles,  V^  V^,  W1 

and  W2  0.0ms 

Particle  diameter  300pm 

Other  Input 

Projectile  frictional  pressure  1M  Pa 

Projectile  mass  2kg 


TABLE  2: 

Input  data  used  in  the  calculations  (source  Refs.  14  and  15 


8. 

MOMENCLATURF 

V 

3  2 

heat  transfer  parameters  for  the  gas  and  particles 

respectively  (see  equations  6.12  and  6.13) 

A 

projectile  cross-sectional  area 

A 

s 

total  surface  area  in  a  control  cell 

At 

Atwood  Number 

b 

burning  rate  of  particles 

C1  ' 

C2 

specific  heats  of  gas  and  particle  respectively 

Cf 

interphase  friction  coefficient 

CF  I  PS 

interphase  momentum  transfer  coefficient 

(see  equation  5.3) 

0 

0 

initial  particle  diameter 

9 

acceleration  due  to  gravity 

h 

enthalpy 

h 

c 

heat  of  combustion 

k 

thermal  conductivity 

*21 

rate  of  interphase  mass  transfer 

n 

propellant  burning  rate  index 

*4  u 

Nusselt  Number 

D 

pressure 

3b 

pressure  at  the  base  of  the  projectile 

3f 

projectile  frictional  pressure 

4 

rate  of  heat  transfer 

radial  co-ordinate 

! 

volume  fraction 

time 

temperature 

ignition 

propellant  ignition  temperature 

particle  surface  temperature 
Velocity  vector 
control  cell  volume 
lateral  co-ordinate 
axial  co-ordinate 

specific  heat  ratio  of  gas 
diffusion  coefficient 
propellant  burning  rate 
proportionality  constant 
wavelength 
dynamic  viscosity 
density 

intergranular  stress 
wall  shear  stress 

pertaining  to  the  gas  phase 
pertaining  to  the  solid  phas 
pertaining  to  the  interface 
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Figure  3  (cont'd)  :  Observations  in  the  experiment  with 

a  tilted  interface  (cont'd) 
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Figure  6  (cont'd) 
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Figure  7  (cont'd)  :  The  effect  of  viscosity 
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Figure  8  (cont'd)  :  The  effect  of  the  wall 

(a)  c  =  1  &  k  =  3 . 9 

(b )  c  =  1 0  &  k  =  3.9 
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Figure  9  (rout'd) 
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Figure  10: 

The  gun  barrel  considered 


Figure  1 1 

The  finite  domain  grid  used 
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The  effect  on  the  pressure 
of  the  number  of  cells  used 
in  the  axial  direction 


Figure  22  The  particle  volume  fraction  in  the  two 

calculations  at  r/R=0.64 

-  without  initial  non-uniformity 

_  with  an  initial  non-uniformity 
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Figure  23 


The  particle  volume  fraction  radial  distribution 
history  at  (a)  Z/L  =  0.33  and  (b)  Z/L  =  0.88 


Figure  24 

Sketch  of  the  growth  of  an  initial  volume 
fraction  non-uniformity  in  a  gun  barrel  situation 


